#include <string>
#include <algorithm>
#include <numeric>
#include <istream>
#include <vector>
#include <fstream>
#include <deque>
#include <map>
#include <stdexcept>
#include <ext/hash_map>


// boost functions
#include <boost/lexical_cast.hpp>

// my functions
#include <portable_timer.h>
#include <progress_indicator.h>
#include <blast_results.h>
#include <print_error.h>
#include <ultoa.h>
#include <index_string.h>

#include "get_arg.h"

#include <gc.h>
using std::cerr;
using std::ifstream;
using std::ofstream;
using std::cout;
using std::string;
using std::vector;
using std::deque;
using std::sort;
using boost::lexical_cast;
using std::ostringstream;





struct  t_gene_pair_blast_results
{
	unsigned prot_id1;
	unsigned prot_id2;
	float	 bitscore;
	float	 e__value;
	unsigned align_len;
	t_gene_pair_blast_results(	unsigned prot_id1_,
								unsigned prot_id2_,
								t_blast_results results):
								prot_id1( prot_id1_),
								prot_id2( prot_id2_),
								bitscore( results.bitscore),
								e__value( results.e__value),
								align_len(results.align_len){}
};

bool operator < (const t_gene_pair_blast_results& a, const t_gene_pair_blast_results& b)
{
	if (a.bitscore < b.bitscore)
		return true;

	// if same evalue
	else if (a.bitscore == a.bitscore)
	{
		//only save if lower e_value
		if (a.e__value > b.e__value)
			return true;
		else if (a.e__value == b.e__value)
		{
			//  or shorter length
			if (a.align_len > b.align_len)
				return true;
		}
	}
	return false;
}

typedef deque<__gnu_cxx::hash_map<unsigned, unsigned> >		t_matrix_blast_result_indices;
typedef __gnu_cxx::hash_map<unsigned, unsigned>::iterator	t_blast_result_index_iterator;
//typedef deque<std::map<unsigned, unsigned> >		t_matrix_blast_result_indices;
//std::map<unsigned, unsigned>::iterator		t_blast_result_index_iterator;

//________________________________________________________________________________________

//		lookup_prot_name

//________________________________________________________________________________________
inline unsigned lookup_prot_name(t_index_string& prot_ids, const string& prot_name)
{
	unsigned prot_id;
	if (!prot_ids.is_indexed(prot_name, prot_id))
		throw std::runtime_error("Protein identifier in blast result [" + prot_name +
								"] not found in the supplied protein id to gene id "
								"lookup table.");
	return prot_id;
}

//________________________________________________________________________________________

//		update_blast_results_matrix

//________________________________________________________________________________________
void update_blast_results_matrix(	deque<t_gene_pair_blast_results>& 	blast_results,
									t_matrix_blast_result_indices&		matrix_blast_results,
									unsigned	  						gene_id1,
									unsigned	  		 			   	gene_id2,
									const t_gene_pair_blast_results&	new_blast_result)
{
	if (matrix_blast_results[gene_id1].count(gene_id2) == 0)
	{
		// save index into new blast results in index
		matrix_blast_results[gene_id1][gene_id2] = blast_results.size();
		blast_results.push_back(new_blast_result);
		assert(blast_results[matrix_blast_results[gene_id1][gene_id2]].prot_id1 ==
					new_blast_result.prot_id1);
		assert(blast_results[matrix_blast_results[gene_id1][gene_id2]].prot_id2 ==
					new_blast_result.prot_id2);
		return;
	}

	unsigned blast_result_index = matrix_blast_results[gene_id1][gene_id2];
	assert(blast_result_index < blast_results.size());
	if (blast_results[blast_result_index] < new_blast_result)
		blast_results[blast_result_index] = new_blast_result;
}


//________________________________________________________________________________________

//		print_saved_homolog_pairs

//________________________________________________________________________________________
void print_saved_homolog_pairs(		const t_index_string&					prot_ids,
									const t_index_string&					gene_ids,
									const deque<t_gene_pair_blast_results>&	blast_results,
									const t_matrix_blast_result_indices& 	matrix_blast_results,
									ofstream&								ofs,
									t_progress_indicator& 					dots)
{
	unsigned cnt_multiple_transcript_genes = 0;
	dots.set_every(100000);
	for (unsigned gene_id1 = 0; gene_id1 < matrix_blast_results.size(); ++gene_id1)
	{
		t_matrix_blast_result_indices::value_type::const_iterator jj = matrix_blast_results[gene_id1].begin();
		t_matrix_blast_result_indices::value_type::const_iterator end = matrix_blast_results[gene_id1].end();
		for (; jj != end; ++jj)
		{
			++dots;
			const t_gene_pair_blast_results& res = blast_results[jj->second];
			unsigned gene_id2 = jj->first;
			ofs	<< gene_ids.string_at(gene_id1) << "\t"
				<< gene_ids.string_at(gene_id2) << "\t"
				<< prot_ids.string_at(res.prot_id1) << "\t"
				<< prot_ids.string_at(res.prot_id2) << "\n";
			++cnt_multiple_transcript_genes;
		}
	}
	assert(cnt_multiple_transcript_genes == blast_results.size());
}







//________________________________________________________________________________________

//		parse_blast_results_from_file

//________________________________________________________________________________________
void parse_blast_results_from_file(	const string&						file_name,
									deque<unsigned>&					prot_id_to_gene_id,
									t_index_string&						prot_ids,
									t_index_string&						gene_ids,

									// multiple transcript genes saved
									deque<t_gene_pair_blast_results>&	blast_results,
									t_matrix_blast_result_indices& 		matrix_blast_results,

									// single transcript genes printed out immediately
									deque<char>&						gene_ids_single_prot,
									unsigned&							cnt_single_transcript_gene_results,
									ofstream&							ofs,
									t_progress_indicator& 				dots)
{
	dots.set_every(25000);

	ifstream is(file_name.c_str());
	if (!is)
		throw std::runtime_error("Could not open " + file_name);

	unsigned cnt_lines = 0;

	try
	{
		// get data for one protein pair
		string prot_name1, prot_name2;
		is >> prot_name1 >> prot_name2;
		t_blast_results new_blast_result;
		is >> new_blast_result;

		while (is.good())
		{
			cnt_lines++;
			dots++;

			// get protein and gene indices
			unsigned prot_id1 = lookup_prot_name(prot_ids, prot_name1);
			unsigned prot_id2 = lookup_prot_name(prot_ids, prot_name2);
			unsigned gene_id1 = prot_id_to_gene_id[prot_id1];
			unsigned gene_id2 = prot_id_to_gene_id[prot_id2];

			// ignore comparisons between transcripts of the same gene
			if (gene_id1 != gene_id2)
			{

				// sort gene indices so that comparisons don't have to be symmetrical
				if (gene_id1 > gene_id2)
				{
					std::swap(gene_id1, gene_id2);
					std::swap(prot_id1, prot_id2);
				}

				// no more transcripts for this gene, just print out
				if (gene_ids_single_prot[gene_id1] && gene_ids_single_prot[gene_id2])
				{
					ofs	<< gene_ids.string_at(gene_id1) << "\t"
						<< gene_ids.string_at(gene_id2) << "\t"
						<< prot_ids.string_at(prot_id1) << "\t"
						<< prot_ids.string_at(prot_id2) << "\n";
					++cnt_single_transcript_gene_results;
				}
				else
				// save if necessary
					update_blast_results_matrix(blast_results,
												matrix_blast_results,
												gene_id1,
												gene_id2,
												t_gene_pair_blast_results(	prot_id1,
																			prot_id2,
																			new_blast_result));
			}


			// get next line
			is >> prot_name1 >> prot_name2;
			is >> new_blast_result;
		}

		if (is.eof())
			return;
	}
	catch(std::runtime_error& e)
	{
		throw std::runtime_error("Invalid blast results format on line #" +
									ultoa(cnt_lines) + " of " + file_name + "\n" + e.what());
	}
	throw std::runtime_error("Invalid blast results format on line #" +
								ultoa(cnt_lines) + " of " + file_name);


}


//________________________________________________________________________________________

//		parse_protid_gene_id_from_file

//________________________________________________________________________________________
void parse_protid_gene_id_from_file(const string&						file_name,
									deque<unsigned>&					prot_id_to_gene_id,
									t_index_string&						prot_ids,
									t_index_string&						gene_ids,
									deque<char>&						gene_ids_single_prot,
									t_progress_indicator& 				dots,
									t_program_args&						args)
{
	ifstream is(file_name.c_str());
	if (!is)
		throw std::runtime_error("Could not open " + file_name);

	unsigned cnt_lines = 0;

	// get data for one protein pair
	string prot_name, gene_name;
	is >> prot_name >> gene_name;
	dots.set_every(5000);
	while (is)
	{
		dots++;
		cnt_lines++;
		bool added;
		unsigned prot_id = prot_ids.index(prot_name, added);
		lookup_prot_name(prot_ids, prot_name);
		// warn duplicatd gene_ids.string_at(prot_id_to_gene_id[prot_id]) and gene_name
		if (added)
		{
			unsigned gene_id = gene_ids.index(gene_name, added);

			// has a transcript for this gene been encountered before.
			// if not, it may be a single transcript gene
			if (added)
			{
				assert(gene_id == gene_ids_single_prot.size());
				gene_ids_single_prot.push_back(1);
			}
			else
				gene_ids_single_prot[gene_id] = 0;

			// save gene_id lookup
			prot_id_to_gene_id.push_back(gene_id);

		}
		else
		{
			args.CERR() << "WARNING:\tDuplicated prot ID " << prot_name
						<< " with previous gene of " << gene_ids.string_at(prot_id_to_gene_id[prot_id])
						<< " and now with gene of " << gene_name << "\n";
		}

		// get next line
		is >> prot_name >> gene_name;
	}

	if (is.eof())
		return;
	throw std::runtime_error("Invalid protein / gene name on line #" +
								ultoa(cnt_lines) + " of " + file_name);


}





//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//		Main

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
int main(int argc, char* argv[])
{
	t_portable_timer timer;

	// don't use C streams
	std::ios::sync_with_stdio(false);

	// saving parameters
	string exe_parameters;
	for (int i = 0; i < argc; ++i)
		exe_parameters += string(argv[i]) + ' ';

	// we use cout to print progress so flush per call
	cout.setf(std::ios_base::unitbuf);

	// get command line arguments
	t_program_args args;
	try
	{
	if (!process_command_line_options(argc, argv, args))
		return 1;

	//
	//	protein to gene id lookup
	//
	deque<unsigned>						prot_id_to_gene_id;
	t_index_string						prot_ids;
	t_index_string						gene_ids;
	deque<char>							gene_ids_single_prot;
	{
		t_progress_indicator 				dots(args.VERBOSE(),  5000);
		parse_protid_gene_id_from_file(	args.prot_to_gene_id,
										prot_id_to_gene_id,
										prot_ids,
										gene_ids,
										gene_ids_single_prot,
										dots,
										args);
		dots.finish();
		args.VERBOSE() << "\t" << gene_ids.size() << " gene identifiers in total\n";
		args.VERBOSE() << "\t" << prot_ids.size() << " protein transcript identifiers in total\n";
		assert(gene_ids.size() == gene_ids_single_prot.size());
		unsigned cnt_single_transcript_gene =
			std::accumulate(gene_ids_single_prot.begin(),  gene_ids_single_prot.end(), 0U);
		args.VERBOSE() << "\t" << cnt_single_transcript_gene << " genes with single transcripts\n";
	}


	//
	// open results file
	//
	ofstream ofs(args.gene_prot_pairs.c_str());
	if (!ofs)
		throw std::runtime_error("Could not open " + args.gene_prot_pairs);

	//
	//	parse blast results
	//
	//
	deque<t_gene_pair_blast_results>	blast_results;
	t_matrix_blast_result_indices 		matrix_blast_results(gene_ids.size());
	unsigned cnt_transcripts = 0;
	unsigned cnt_single_transcript_gene_results;
	{
		args.VERBOSE() << "\tParsing blast results from " << args.input_files.size() << " files\n";
		t_progress_indicator 				dots(args.VERBOSE(),  5000);
		dots.use_timer();

		for (unsigned i = 0; i < args.input_files.size(); ++i)
		{
			parse_blast_results_from_file(	args.input_files[i],
											prot_id_to_gene_id,
											prot_ids,
											gene_ids,
											blast_results,
											matrix_blast_results,
											gene_ids_single_prot,
											cnt_single_transcript_gene_results,
											ofs,
											dots);
			if (i % 10 == 9)
				dots.print_indicator(ultoa(i+1));
			//args.VERBOSE() << "\n" << dots.value() << "," << blast_results.size() << "\n";
		}

		cnt_transcripts = dots.value();
	}







	// print out rest of results (those from genes with multiple transcripts
	{
		args.VERBOSE() << "\tPrinting out gene/protein homologue pairs for genes with multiple transcripts\n";
		t_progress_indicator dots(args.VERBOSE(),  5000);
		print_saved_homolog_pairs(	prot_ids,
									gene_ids,
									blast_results,
									matrix_blast_results,
									ofs,
									dots);
	}

	args.VERBOSE() << "\t" << string(60, '_') << "\n";
	args.VERBOSE() << "\t" << gene_ids.size() << " genes in total\n";
	args.VERBOSE() << "\t" << cnt_single_transcript_gene_results << " blast results from genes with single transcripts\n";
	args.VERBOSE() << "\t" << blast_results.size() << " blast results for multiple transcript genes\n";
	args.VERBOSE() << "\t" << timer.reset() << " Seconds for " << cnt_transcripts << " parsing blast comparisons\n";
	}
	catch (std::exception& e)
	{
		std_print_error(args.CERR(), e.what());
		return 1;
	}

	return 0;
}



